Contact Pair Dynamics During Folding of a Model Globular Protein, Hp-36 
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The dynamics of contact pair formation between various hydrophobic residues during folding of 
a model protein Hp-36 is investigated by Brownian dynamics simulation. Hydropathy scale and 
non-local helix propensity of amino acids are used to model the complex interaction potential. The 
resulting structure of the model protein mimics the native state of the real protein with a RMSD of 
4.5 A. A contact pair distance time correlation function (CPCF), Cp(t), is introduced which shows 
multistage decay, including a slow late stage dynamics for a few specific pairs. These pairs determine 
the long time folding rate. Dynamics can be correlated with the landscape, relative contact order 
and topological contact. 

PACS numbers: 87.15.Aa,87.15.Cc,87.15.He,87.15.-v,83.10.Mj 



The dynamics of folding of an extended protein chain 
at high temperature (or high urea concentration) to its 
unique folded state at low temperature (or low urea con- 
centration) is a highly complex problem with many in- 
teresting aspects. Recent experimental, theoretical and 
computer simulations studies [1 — 7] have unearthed and 
explained many fascinating aspects of folding, although 
still many others remain to be explored. The paradigm of 
landscape (with the idea of folding funnel) has provided 
new insight into the problem |^|, ^). Experimental data 
on the rate of folding of a large number of small proteins 
have suggested a close relation between the the relative 
contact order and the rate of folding 0. The relative 
contact order denotes the average sequence distance be- 
tween the hydrophobic pair contacts and is defined as M , 
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where (i,j) are the specific hydrophobic pair contacts, N c 
is the number of contacts while L is total number of hy- 
drophobic amino acids present in the protein. Sj and 
are the sequence number along the contour of the chain. 
The rate of folding was found to decrease nearly expo- 
nentially with RCO. The dynamics of such non-local 
contact formation holds the key to the understanding of 
dynamics of folding. However, this aspect has remained 
largely unexplored. 

An attractive way to explore pair dynamics of hy- 
drophobic contact is via the technique of fluorescence 
resonance energy transfer (FRET), fn FRET, one mea- 
sures the time dependence of energy transfer from a cho- 
sen donor fluorophore to a chosen acceptor. The rate of 
transfer may be due to dipolar interactions and the rate 
of transfer is given by the well known Forster expression 
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where k ra( i is the radiative rate and Rp is the Forster 
radius. By suitably choosing donor-acceptor pair, Rp 
can be varied over a wide range. This allows the study of 
the dynamics of pair separation, essential to understand 
protein folding JlO| . k ra d is typically less than (but of the 
order of) 10 9 sec -1 . Thus Forster transfer provides us 
with a sufficiently fast camera to take snapshots of the 
dynamics of contact pair formation. 

In this study, we have studied contact formation by 
Brownian dynamics simulations of a model protein Hp- 
36 which is one of the smallest protein that folds au- 
tonomously to a stable compact ordered structure, with a 
large helix content jll) . Hp-36 is a subdomain of chicken 
villin which is implicated in the formation of microvilli 
in the absorbtive epithelium of the gut and the proximal 
tube of the kidney All atom simulation study on 

this Hp-36 have revealed at least two pathways of folding 
p3[ . Earlier, several studies of Hp-36 were presented us- 
ing Monte Carlo technique |l4|] and Brownian dynamics 

The model studied here is constructed by taking two 
atoms for a particular amino acid. The smaller atom rep- 
resents the backbone C a atom of real protein while the 
bigger atom mimics the whole side chain residue. Con- 
struction of the model protein has been described in de- 
tail elsewhere Similar types of model (with more 
rigorous force field) have been introduced by Scheraga et 
al. recently fl7j] . The total potential energy function of 
the model protein Vrotal is written as, 



V T otal = V B +Ve+V T + VLJ + Vheli. 
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where Vb and V$ are the potential contributions due to 
vibration of bonds and bending motions of the bond 
angles. Standard harmonic potential is assumed for 
the above two potentials with spring constants 43.0 kJ 
mole _1 A -2 and 8.6 kJ mole _1 A -2 for the bonds between 
backbone atoms and bonds joining side residues with the 
backbone atoms, respectively. In case of the bending 
potential, spring constant is taken to be 10.0 kJ mol -1 
rad -2 . V T (= £t E^l 1 / 2 )!! + cos(3aV)]) is taken as tor- 
sional potentials for the rotations of the bonds, ct = 1 
kJ mol -1 . The nonbonding potential Vlj is the sum of 
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FIG. 1: (a) The backbone structure of the model protein with 
the lowest RMSD (4.5 A), (b) The backbone structure of the 
native state of real Hp-36. 



the pair interactions between the atoms and is given by, 



V LJ = 4j2ti 



r / \ 12 













(4) 



where r,j and 



are the distance and interaction be- 



tween the i-th and j-th atom. Oi 



j) and 



Sizes and interactions are taken to be the 



same (1.8 A and 0.05 kJ mol -1 , respectively) for all the 
backbone atoms as they represent the C a atoms in case of 
real proteins. Side residues, on the other hand, carry the 
characteristics of a particular amino acid. Different sizes 
of the side residues are taken from the values given by 
Levitt Q . Interactions of the side residues are obtained 
from the hydrophobicities of the amino acids. We con- 
struct effective potential guided by the well-known statis- 
tical mechanical relation between potential of mean force 



and radial distribution function as Vi 
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JL9| . Strong correlation among the hydrophobic groups 
(absent among the hydrophilic amino acids) implies that 
the hydrophobic amino acids should have stronger effec- 
tive interaction than the hydrophilic groups. So the in- 
teraction parameters of the side residues can be mapped 
from the hydropathy scale p0| by using a linear equation 
as given below, 
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where, tu is the interaction parameter of the zth 
amino acid with itself. e m i n (=0.2 kJ mol -1 ) and 
£max( = H-0 kJ mol -1 ) are the minimum and maximum 
value of the interaction strength chosen for the most 



hydrophilic(arginine) and most hydrophobic (isoleucine) 
amino acids, respectively. Ha is the hydropathy in- 
dex of zth amino acid given by Kyte and Doolittle p(J 
and H m i n (=-4.5) and H max {=4.5) are the minimum and 
maximum hydropathy index among all the amino acids. 
Further details are available in Ref. 16. An important 
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FIG. 2: (a) The multistage temporal evolution of energy En 
and minimized energy E^ rn (left ordinate) , radius of gyration 
R-y (right ordinate) and RMSD (right ordinate in parenthe- 
sis) are plotted in the same time axis to show the dynamical 
correlation, (b) Dynamics of the topological contact Nt opo 
(left ordinate) and the RCO (right ordinate) are plotted in 
the same time axis which show a similar dynamical growth. 

part of secondary structure of the real protein is the for- 
mation of a helix. In the absence of hydrogen bonding, 
we introduce the following effective potential among the 
backbone atoms to mimic the helix formation along the 
chain of residues that show high helix propensity, 
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where r^t+a and r^f+a are the distances of iih atom 
with i + 2 and i + 3 th atoms, respectively, r/, is 
the equilibrium distance and is taken as 5.5 A, moti- 
vated by the observation that the distance of with 
ri + 2 and ri+$ is nearly constant at 5.5 A in an a he- 
lix. The summation excludes the first and last three 
amino acids as there is less helix formation observed in 
the ends of the protein chain pi]] . The force constant 
for the above harmonic potential is mapped from the 
helix propensities Hpi taken from Scholtz et al. [ p2[ , 

^■i ^-alanine Hp^X (/C alanine ^-glycine) • ^alanine and 

^■glycine are the force constants for alanine and glycine, 
17.2 and 0.0 kJ mol -1 , respectively. Next, the influ- 
ence of the neighboring amino acids for the formation 
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of helix has been considered by taking an average of the 
spring constants as K\~ z = ^[JCi + + JCt+2] and 
— j[JCi + JCi + i + /Cj + 2 + /Q+3], with the condi- 
tion that ifj 1-3 , > as the force constant must 
remain positive. The above formulation of helix poten- 
tial is motivated by the work of Chou and Fasman about 
the prediction of helix formation that the neighbors of a 
particular amino acid should be considered rather than 
its own helix propensity [p3| . 

Figure 1(a) shows the best folded structure obtained in 
our simulations while figure 1(6) shows the real one. The 
RMSD calculated over backbone with the real native 
Hp-36 is 4.5 A which is reasonable for such a simplified 
model. The initial configuration of the model protein 
was generated by configurational bias Monte Carlo tech- 
nique J24|] . Atoms attached to a single branch point were 
generated simultaneously. Then the initial configuration 
was subjected to Brownian dynamics simulation for the 
study the folding. Time evolution of the model protein 
was carried out according to the motion of each bead as 
below, 

Ti(t + At) = Ti (t) + ^F<(t)Ai + Arf (7) 

where each component of Ar^ is taken from a Gaussian 
distribution with mean zero and variance ({Arf a ) 2 ) — 
2D At [2^]. Ti(t) is the position of the ith atom at 
time t and the systematic force on ith atom at time t is 
Fj(i). The time step At is taken as 0.001. D% is the dif- 
fusion coefficient of the i-th particle calculated from the 
Stokes- Einstein relation Z)j = ^ jB [. . Ri is the radius of 
the i-th atom and rj is the viscosity of the solvent, kp and 
T are the Boltzmann constant and temperature, respec- 
tively. Simulations have been carried out for Af number 
of different initial configurations, where AT = 584. 

The potential energy En, the radius of gyration i? 7 
and the RMSD (calculated over backbone from the real 
native structure of Hp-36 obtained from protein data 
bank |ll], ^6)) all exhibit the multistage dynamics as 
shown in figure 2(a) in the same time axis. There is an 
initial sharp hydrophobic collapse followed by a slower 
decay is observed till 500r (w 200 ns). A long plateau 
follows in the final stage which exists for a very long 
time ( 2000r w 1/is) before the protein reaches its final 
lowest energy state at around 2400r. Energy values of 
the corresponding inherent structures obtained by con- 
jugate gradient is shown by the symbols which depicts 
the decreasing local minima attained by the system un- 
til it reaches the final folded state. Total number of hy- 
drophobic topological contact N t0 p and relative contact 
order RCO are plotted in figure 2(b). N top o is defined to 
be formed if two hydrophobic side chain residues come 
within a distance of 8.5 A. RCO is calculated from Eq. 
1. Both N topo and RCO show a multistage increase sig- 
nifying the participation of nonlocal contacts. 

Folding can be probed microscopically by monitoring 
the dynamics of separation between different amino acid 
pairs. The widely different time scales of movement of 
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FIG. 3: Dynamics of contact formation of different side 
residues with the 9-th side residue is shown. The multistage 
relaxation process in the dynamical quantities originates from 
the diverse dynamics of the contact pairs. 



all the different pairs together give rise to an overall pic- 
ture of the dynamics of folding which is reflected in the 
macroscopic quantities. The effective dynamics of pair 
separation can be described by introducing a new pair 
correlation function defined below Jlql, 



d ij {t) -tfc'(oo) 
d^'(0) - d^(oo) 
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where, d u (t) = |r.;(t) — r 3 -(i)|. and r,- are the positions 
of the z-th and j-th atom, respectively. Figure 3 shows 
the Cp(t) of the 9th side residue with many other hy- 
drophobic side residues. The three stages of the folding 
process are reflected by mainly three different dynamical 
behavior seen amongst the side residues. Side residues 
closed to the tagged one collapse very fast. Some show 
an initial shoulder and only a few show the plateau in the 
long time decay that correlates with the similar plateau 
observed in case of other dynamical quantities. The fi- 
nal decrease in energy is observed when all the different 
contact pair correlation functions decay at around 2400r. 

We have calculated survival probability Sp(t) in FRET 
using Forster energy transfer rate from Eq. 2 for the 
9 — 35 pair along the Brownian dynamics trajectory lead- 
ing to the most stable structure. The Forster radius Rf 
is taken as 10 A. Figure 4 shows initial very slow decrease 
of Sp(t) to be followed by a sudden drop at around 2400t 
as seen in case of the different dynamical properties dis- 
cussed above. Sp(t) is found to be relatively insensitive 



to k r 



The resemblance with the real native state of the 



protein and the dynamics of RMSD show the validity of 
the model used for this protein. Contact pair dynamics 
and the time evolution of energy, radius of gyration, rela- 
tive contact order formation etc. brings out the rich and 
diverse dynamics of protein folding. The initial ultrafast 
hydrophobic collapse signify that the upper part of the 
funnel is steep - followed by a change in slope. Rate de- 
termining step, however, arises from the final stage of 
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folding on a very flat and rugged underlying landscape 
marked by the large conformational entropy barrier with 
little energy change This entropic bottleneck arises 
from the necessity to form long range hydrophobic con- 
tacts, as envisaged by Dill and Wolynes. The atoms mim- 
icking the whole side chain of the real protein play a very 
important role for structural and dynamical aspects in 
this study. Moreover, the new contact pair correlation 
function and FRET probes the folding events in minute 
detail. 
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FIG. 4: FRET survival probability for R F = 10.0 A is plotted 
for different radiative rates to capture the wide variety of time 
scales involved in FRET. 
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